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I.   INTRODUCTION 

For  a  structural  analysis  of  any  non-trivial  elasto-static 
problem  by  the  finite  element  method,  a  large  amount  of  data  must 
be  provided  to  the  computer  system  which  performs  the  calculations. 
The  required  information  includes  the  arrangement  of  nodes  within 
each  element,  the  coordinates  of  all  nodes,  and  the  applied  loads. 
After  the  numerical  values  of  this  data  have  been  determined,  which 
is  a  large  problem  in  itself,  the  information  must  be  supplied  to 
the  computer  system.   This  is  normally  accomplished  by  punching  the 
appropriate  data  on  computer  cards.   The  typing  of  these  cards  is 
also  a  difficult  and  time-consuming  job.   More  important,  however, 
is  the  ever-present  possibility  of  human  errors  during  the  performance 
of  these  tasks.   If  an  error  is  not  detected  prior  to  submitting 
the  data  for  processing,  expensive  computer  time  would  be  wasted. 

With  the  above  considerations  in  mind,  it  is  readily  apparent 
that  a  computer  system,  which  would  generate  most  of  the  data  required 
by  the  primary  computer  system  for  a  complete  structural  analysis, 
would  be  very  useful.   In  addition,  this  mesh  generator  should 
require  as  little  input  data  as  possible. 

References  5,  8,  and  10  contain  reports  of  efforts  which  have 
been  directed  at  coding  computer  systems  which  will  aid  the  analyst 
in  the  interpretation  of  output  data.   A  few  computer  systems,  which 
are  reported  in  Refs.  1,  2,  and  9,  have  been  coded  to  automatically 
generate  input  data.   However,  these  input  data  generation  systems 


are  either  for  types  of  finite  elements  other  than  isoparametric 
elements  or  for  two-dimensional  elements. 

Two  input  data  generator  systems  for  three-dimensional  iso- 
parametric finite  elements  will  be  presented  in  this  thesis.   One 
system  is  for  quadratic  elements  and  the:  second  system  is  for  cubic 
elements . 


II.   ISOPARAMETRIC  FINITE  ELEMENTS 

The  isoparametric  concept  will  be  discussed  in  this  chapter  in 
sufficient  detail  to  establish  only  the  procedures  for  constructing 
isoparametric  elements  and  the  validity  of  relationships  used  to 
obtain  the  consistent  gravity  and  pressure  loads.   The  readers  are 
referred  lo  Refs.  4  and  7  for  further  details  of  the  isoparametric 
concept . 

A.   DISPLACEMENT  FUNCTIONS 

The  displacements  of  nodes  in  any  one  element  are  obtained  by 

using  shape  functions  which  are  defined  in  a  system  of  dimensionless 

coordinates  (£>,T|,0  that  are  unique  for  that  element.   The  coordinates 

range  from  minus  one  to  plus  one. 

u(x,  y,  z)  =  SN.(§,71,0  u. 

v(x,  y,  z)  =  EN .(§,T|,0  v 

111  (1) 

w(x,  y,  z)  =  £N.(§,T1,0  w. 

i  *■  ± 

where  N.  (^j^jQ  are  the  shape  functions  and  u.,  v.,  w.  are  the  dis- 
placement components  of  node  "i"  in  a  global  reference  system  (x ,  y,  z) , 
where  "i"  assumes  values  from  1  to  the  total  number  of  nodes  in  the 
element.   This  global  reference  system  is  then  related  to  the  dimen- 
sionless coordinates  by  the  same  shape  functions  as  follows: 
x(S,T],0  =  SN.(S,T|,0  Xi 

y(S,T|,Q  =  Sn .<M.C)  y± 

11  (2) 

z(§,^i,o  =  sn  a/i.o  z. 

i   i  l 
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where  x.,  v.,  z.  are  the  coordinates  of  node  "i"  in  the  global 
reference  system  [4], 

The  shape  functions  referred  to  above  have  been  obtained  from 
[73.   They  are  listed  in  section  2.C. 


B.   TRANSFORMATIONS 

In  order  to  obtain  the  consistent  gravity  and  pressure  loads, 
the  volumetric  increment,  dxdydz ,  must  be  found  in  terms  of  the 
dimensionless  coordinates.   That  operation  requires  the  Jacobian 
matrix  [4], 

ox/S§  By/3§  Bz/S§ 

[J]  -      ox/o7]   dy/aTl  dz/57] 

ox/dC  oy/dC  oz/BC 


(3) 


The  element  of  volume  is  then  transformed  as  follows 
dV  =  det[j]d£  dT]  dC 


(4) 


C.   SHAPE  FUNCTIONS 

In  this  listing  of  the  shape  functions,  attention  is  directed 
to  the  tri-quadra tic  and  tri-cubic  elements,  for  which  the  mesh 
generator  programs  were  written. 

1.   Quadratic  Elements 

See  Figure  1  for  identification  of  the  nodes. 
Corner  nodes:   1,  3,  5,  l\    13,  15,  17,  and  19. 


N.  =  (1/8)(1  +  §  )(1  +T\   )(l  +  C  )(§   +  T]  +  C   -  2) 

1  o       oy  o    o    o    o 

where   %,      =  §.§   and  §.    is   _    1;    similarly    for   T|     and   C,    . 

Mid-side  nodes:      2,    6,    14,    and    18. 

N.    =    (1/4)  (I   -   §2)(1  +  71   )(1   +  C   ) 
1  00 


(5) 


(6) 
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Mid-side  nodes:   4,  8,  16,  and  20. 

N.  =  (1/4)  (1  -  II2)  (I  +  §o)(l  +  CQ)  (7) 

Mid-side  nodes:   9,  10,  11,  and  12.' 

N.  =  (1/4)  (1  -  C2)(l  +  S0)(l  +  \)  (8) 

2 .   Cubic  Elements 

See  Figure  2  for  identification  of  nodes. 

Corner  nodes:   1,  4,  7,  10,  21,  24,  27,  and  30. 

N.  =  (1/64)  (1  +  §o)(l  +  Tlo)(l  +  C0)[9(§2  +T\2  +  C2)    -  19]  (9) 

Mid-side  nodes:   2,  3,  8,  9,  22,  23,  28,  and  29. 

N .  =  (9/64)  (1  -  §2)  (1  +  9§o)  (1  +  T]o)  (1  +  C0)  (10) 

where  § .  =  -  1/3 ,  T] .  =  +  1 ,  and  C.  =  +  1 . 
1  x  1 

Mid-side  nodes:   5,  6,  11,  12,  25,  26,  31,  and  32. 

N.  =  (9/64)  (1  -  T]2)(l  +  971  )(1  +  §  )  (1  +  C  )  (11) 

i  o       o       o 

where  T\.    =  t   1/3,  §.  =  ±  1,  C-  -  1*  1- 

Mid-side  nodes:   13,  14,  15,  16,  17,  18,  19,  and  20. 

N.  =  (9/64)  (1  -  C2)(l  +  %   )(1  +  5  )(1  +  71  )  (12) 

1  o        o        o 


where  £,  =  -  1/3,  §.  =  J"  1,  and  71.  -  -  1. 


Note  that  the  convention  for  numbering  of  nodes  within  any 
element  is  to  begin  at  the  location  where  §  =  7|  =  £  =  1  and  to 
consecutively  number  the  nodes,  proceeding  in  a  counter-clockwise 
direction  around  the  C-axis  (see  Figures  1  and  2). 
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Figure  1.   Quadratic  Element 


13 


Figure  2.   Cubic  Element 
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III.   DISCUSSION  OF  COMPUTER  PROGRAMS 

The  computer  programs  presented  in  this  thesis  were  written  to 
support  a  computer  system  called  TRISOP.   TRISOP  was  coded  by  Professor 
G.  Cantin  of  the  U.  S.  Naval  Postgraduate  School,  Monterey,  California. 
TRISOP  performs  a  structural  analysis  of  three-dimensional  elasto- 
static  problems  using  isoparametric  finite  elements.   Information 
required  by  TRISOP  for  any  problem  includes  general  mesh  parameters, 
title  cards,  format  statements,  and  material  properties;   Young's 
modulus  of  elasticity,  Poisson's  ratio,  and  the  coefficient  of  thermal 
expansion.   Fifteen  cards  are  needed  to  supply  this  data  to  TRISOP. 
The  bulk  of  the  data  required  by  TRISOP  is  the  element  connectivity, 
structure  geometry,  loading,  and  boundary  conditions.   The  boundary 
conditions  will  vary  with  each  problem;  therefore,  no  effort  has 
been  made  to  automate  this  process. 

The  mesh  generating  computer  programs  discussed  in  this  chapter 
compute  the  remainder  of  the  data  required  by  TRISOP:   element 
connectivity,  coordinates  of  nodal  points,  consistent  gravity  load 
vector,  and  consistent  pressure  load  matrix.   Additional  information 
which  is  helpful  in  preparing  an  input  deck  of  cards  for  TRISOP  is  . 
also  computed.   One  of  these  programs  generates  data  for  quadratic 
elements  and  the  second  program  generates  data  for  cubic  elements. 
These  programs  will  be  called  QUAMEG  and  CUMEG. 

A.   ELEMENT  CONNECTIVITY 

A  structure  of  any  shape  can  be  visualized  as  a  cubical  piece 
of  plastic  material  which  can  be  stretched,  compressed,  and  otherwise 
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appropriately  deformed  in  such  a  way  that  it  will  assume  the  geometry 
of  the  structure.   The  element  connectivity  is  determined  when  the 
structure  is  in  this  initial  abstract  cubical  shape.   The  dimensionless 
coordinates  of  this  abstract  structure  are  §  ,  T)  ,  C,    . 

Numbering  of  nodes  is  begun  in  a  face  which  will  have  the 
fewest  nodes.   In  that  face,  numbering  proceeds  sequentially  along 
the  edge  having  the  smallest  number  of  nodes  as  shown  in  Figure  3. 
This  convention  reduces  the  size  of  the  half-bandwidth  of  the  system, 
which  directly  affects  the  space  and  time  required  by  the  computer 
to  solve  the  problem.   The  coordinate  system  is  positioned  such 
that  node  one  is  at  §  =  0,  T|   =1,  and  Q     =1.   An  example  is 
discussed  in  Appendix  A. 

For  this  mesh,  the  connectivities  of  the  elements  follow  simple 
but  tedious  rules  that  were  coded  directly. 

B.   SOLUTION  TIME  AND  SPACE  REQUIREMENTS 

The  time  required  to  solve  the  system  of  equations  in  the  problem 
and  the  space  required  by  the  computer  to  store  all  data  are  important 
information.   They  must  be  determined  before  the  problem  is  submitted 
for  processing.   The  calculation  of  this  information  requires  no 
additional  input  data  and  it  is  therefore  included  in  QUAMEG  and  CUMEG 
to  reduce  the  effort  required  of  the  programmer. 

The  time  required  to  solve  the  equations  is  computed  using 
an  empirical  formula  which  was  constructed  by  Cantin  [4],  This 
calculation  has  been  included  in  the  coding  of  QUAMEG  and  CUMEG. 

The  space  required  to  store  the  various  bits  of  information 
for  a  problem  depends  upon  the  total  number  of  records  and  the 
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length  of  each  record.   These  will  vary  with  different  types  of 
computers.   Simple  formulas  which  generate  the  correct  results  for 
the  IBM  360/67  computer  used  at  the  U.  S.  Naval  Postgraduate  School 
are  employed  in  the  coding  of  QUAMEG  and  CUMEG. 

C.   NODES  ON  BOUNDARY  SEGMENTS 

Before  the  coordinates  of  nodes  which  are  within  the  boundaries 
of  a  structure  can  be  computed,  the  coordinates  of  all  corner  nodes 
on  the  edges  of  the  bounding  surfaces  must  be  specified.   This  problem 
is  simplified  by  the  method  discussed  below. 

The  term  "boundary  segment"  must  first  be  defined  for  the  purpose 
of  this  discussion.   A  segment  must  satisfy  a  few  simple  conditions 
which  depend  upon  the  orientation  of  the  proposed  segment  in  the 
structure  3rm  upon  tuG  type  of  cooruj_nate  system  useu  to  utsciiue 
the  geometry  of  the  structure,.   QUAMEG  and  CUMEG  are  coded  to  allow 
the  use  of  either  rectangular  or  cylindrical  coordinate  systems.   A 
combination  of  the  two  coordinate  systems  in  any  one  problem  is  not 
permitted . 

1.  Rectangular  Coordinates 

When  rectangular  coordinates  are  used,  a  segment  is  any 
straight  line  between  two  corner  nodes  of  the  mesh.   The  element 
divisions  along  this  line  must  be  uniform.   There  are  no  other  restric- 
tions imposed  upon  a  segment  when  rectangular  coordinates  are  used. 

2.  Cylindrical  Coordinates 

When  cylindrical  coordinates  are  used,  the  radius  and  reference 
angle  are  taken  in  planes  parallel  to  the  x-y  plane.   In  this  case 
a  segment  can  be  an  arc  of  circle  centered  on  z  and  in  a  plane 


parallel  to  plane  x-y,  or  a  straight  line  in  a  plane  containing 
the  z-axis.  Again,  element  divisions  along  these  lines  must  be 
uniform. 

In  order  to  determine  the  coordinates  of  nodes  along  a  boundary 
segment,  the  coordinates  of  the  corner  nodes  at  the  ends  of  the 
segment  must  be  provided  in  the  input  data.   In  addition,  the  location 
of  the  segment  within  the  mesh  must  be  specified.   The  method  of 
specifying  the  location  of  a  segment  is  described  in  sections  4.B.3 
and  4.B.4. 

In  Figure  4,  the  vector  AB  is  easily  obtained  from  the  coordinates 
of  nodes  A  and  B.   After  calculating  the  values  of  increments  of 
distances  along  x,  y,  and  z,  the  intermediate  nodes  are  easily 
obtained. 

If  the  bounding  surfaces  of  a  structure  are  of  an  irregular 
configuration  such  that  portions  of  the  edges  of  that  surface  cannot 
be  specified  by  segments,  the  coordinates  of  all  nodes  along  this 
irregular  portion  of  the  bounding  surface  must  be  provided  in  the 
input  data . 

D.   CORNER  NODES 

Once  the  coordinates  of  the  nodes  on  all  of  the  boundary  segments 
are  established,  the  coordinates  of  corner  nodes  on  the  bounding 
surfaces  and  on  interior  slices  are  calculated  using  a  method  described 
below.   The  bounding  surfaces  are  considered  in  the  following  order 
(Figure  3):   face  1  (C'=l)  ,  face  2  (C'=0),  face  3  (§'=0),  face  4  (§'=1), 
face  5  (T)  =1)  ,  face  6  (Tj  =0)  .   Slices  which  are  each  of  constant  Q 
are  then  considered  until  the  coordinates  of  all  corner  nodes  are 
calculated. 
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Figure  4.   Boundary  Segment 
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The  method  used  for  calculating  the  coordinates  of  corner  nodes 
within  the  edges  of  a  bounding  surface  or  slice  can  best  be  understood 
with  the  example  illustrated  in  Figure  5.   In  this  example,  the 
surface  ABCD  is  the  real  surface  for  which  a  mesh  is  desired.   The 
mesh  is  first  constructed  in  the  a$    plane  as  shown. 

At  this  time,  the  coordinates  of  all  the  corner  nodes  on  the 
boundary  are  known.  If  the  coordinates  of  the  interior  nodes  are 
required  to  satisfy  the  Laplacian  equation: 

V2(x)  =  0 

V2(y)  =  0 

V2(z)  =  0  (l3> 

2    2    2    2   2 
where  V   =  0  /£«  +  3  /0£ 

a  set  of  acceptable  interior  nodes  will  be  obtained.   The  dummy 

variables  (#,8)  are  replaced  by  appropriate  combinations  of  §  ,  T|  , 

and  £  ,  depending  upon  the  location  of  the  face  being  considered. 

For  this  computation,  the  standard  Laplacian  molecule  of  finite 

difference  theory  is  used  and  the  final  coordinates  are  obtained 

by  relaxation.   Since  the  coordinates  of  the  nodes  within  the  edges 

of  a  face  are  arbitrary,  the  relaxation  process  can  be  stopped  at 

any  time,  depending  upon  the  degree  of  approximation  desired.   In 

this  code,  sixty  iterations  are  used. 

E.   MID- SIDE  NODES 

After  the  coordinates  of  all  corner  nodes  have  been  established, 
the  coordinates  of  the  mid-side  nodes  are  computed  by  averaging 
the  coordinates  of  adjacent  corner  nodes.   Although  the  mid-side 
nodes  need  not  be  located  exactly  mid-way  between  the  corner  nodes, 
it  is  the  simplest  means  of  determining  their  coordinates. 
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Figure  5 
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F.   CONSISTENT  GRAVITY  LOAD  VECTOR 

The  z-axis  has  been  chosen  as  the  line  of  action  o£  the  gravity 
load.   The  work  done  by  the  force  due  to  gravitational  acceleration 

on  a  body  is 

o  (14) 

U  =  \v  pgwdV 

.here  D  is  the  work  done;  p  is  the  .ass  density  of  the  hody;  g  is 
gravitational  acceleration-,  .  is  the  displaced  of  the  body  in 
the  Z-direction;  and  dV  is  the  element  of  volume. 

*.„a    fnrrpq  actine  at  each  node  in 
The  work  done  by  concentrated  rorces  acting 

the  structure  mesh  is,  in  matrix  notation: 

f  \  (l5) 

vhere  <V>  is  the  consistent  gravity  load  vector  and  ^  are  the 
respective  node  displacements. 

From  Chapter  2,  recall  the  following  two  expressions: 

v(x,y,z)  -p.CS/TU)^ 

dV  =  det[j]  d^dT^dC. 
When  these  expressions  are  substituted  into  equation  14  and  the 
result  identified  to  equation  15,  the  consistent  gravity  load  vector 

is  obtained. 

<V>-pg\U  det[J](N.)d?dW 

-1-1 -I 
In  this  form,  equation  16  is  ideally  suited  for  solution  by 
Guassian  cubature.   In  QUAlffiG,  four  Gauss  points  are  used  in  each 
o£  the  three  directions,  for  a  total  of  64  points  within  each  element. 
Five  Gauss  points  are  used  in  GUlffiG .   The  weight  factors  and  coor- 
dinate values  are  obtained  from  Ref.  11.   The  Jacobian  matrix  is 
determined  by  using  three  function  statements  and  the  shape  functions 

23 


are  determined  with  two  function  statements  in  Subroutine  GRAP 
(Appendix  C)  .   These  function  statements  are  permuted  to  obtain  the 
partial  derivatives  of  the  shape  functions  and  the  values  of  the 
shape  functions  for  all  nodes  in  the  mesh. 

G.   CONSISTENT  PRESSURE  LOAD  MATRIX 

The  consistent  pressure  load  matrix  for  a  uniformly  distributed 
pressure  acting  on  any  face  of  an  element  is  determined  in  a  manner 
similar  to  that  used  for  the  gravity  load. 

The  work  done  by  the  uniformly  distributed  pressure  is: 

U  =  ^A(~pn  dA).(u)  (17) 

where  n  is  a  unit  normal  vector  (positive  outward)  and  "u  is  the 
displacement  vector  at  the  surface  (see  Figure  6).   The  displacement 
vector  u  can  also  be  written  as: 

u  =  u(>:,y,z)T  +  v(x,y,z)j  +  w(x,y,z)k  (18) 

The  vector  ndA  is  found  from  differential  geometry: 
ndA  =  (37/otv)X(d7/3e)  dadB  (19) 

where  r"  =  xi  +  yj  +  zk  is  a  position  vector  of  a  point  on  the  pres- 
surized surface  and  a   and  8  are  dummy  variables,  standing  for  either 

§,  T|,  or  C. 

The  terms  drVcKy  and  dr"/d8  are  simply  two  rows -of  the  Jacobian 
matrix  evaluated  on  the  surface.   After  making  appropriate  substi- 
tutions and  performing  a  few  matrix  manipulations,  the  expression 
reduces  to: 

U  =  \j  ^  <f(cv,8)>  <u*>  dadB  (20) 

-1-1        L  LJ 


where  su.>  is  the  vector  of  nodal  coordinates  for  the  eight  nodes 
on  the  face  and  f(cy,6)  is  a  row  vector  of  functions  of  Q>   and  B. 
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The  integration  is  performed  by  Gaussian  quadrature  and  the 
resultant  vector  is  exhibited  in  convenient  matrix  form.   Five  Gauss 
points  are  used  in  both  QUAMEG  and  CUMEG. 

The  Jacobian  matrix  and  shape  functions  are  determined  in 
Subroutine  GRAP  by  using  the  same  function  statements  as  are  used 
in  the  calculation  of  the  consistent  gravity  load  vector.   The  face 
number  of  the  element  on  which  the  pressure  is  acting  must  be  included 
in  the  input  data.   With  this  information  available,  a  "computed 
go  to"  statement  is  used  in  the  code  to  determine  which  variables 
(5>  TL  C)  will  replace  a  and  (3.   Referring  to  -Figure  6,  the  faces 
are  located  as  described  below. 

Face   Location   a  j3       Face   Location   0/   fi 

1  §  -  +1     Tl  C 

2  Tl  =  +1     C   § 

3  C  -  +1     ?  Tl 

H.   PLOT  OF  STRUCTURE  MESH 

Once  the  coordinates  of  all  nodes  in  the  mesh  have  been  established, 
the  mesh  must  be  inspected  to  determine  if  it  is  acceptable.   Without 
a  plot  of  the  coordinates,  this  task  is  nearly  impossible.   Plotting 
the  mesh  by  hand  is  prohibitively  time-consuming.   Subroutine  GRID 
has  been  included  in  QUAMEG  and  CUMEG  to  generate  an  off-line  printer 
plot  of  the  mesh.   The  parameters  required  for  IBM  360  source  library 
program  DRAW  are  determined  in  GRID.   DRAW  is  called  to  perform  the 
actual  plotting  of  the  mesh.   GRID  generates  a  two-dimensional 
plot  of  the  mesh  on  slices  at  constant  values  of  £  .   For  the  example 
shown  in  Figure  3,  faces  1  and  2  and  five  slices  at  constant  values 
of  C,      would  be  plotted. 
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The  user  can  specify  that  the  scaling  be  performed  automatically 
by  DRAW.   However,  that  choice  could  result  in  the  x-scale  being  not 
equal  to  the  y-scale.   The  plot  of  the  mesh  would  be  distorted  in 
that  case.   The  scaling  is  done  in  GRID  to  insure  that  the  scales 
will  be  equal  and  that  the  smallest  possible  scale  is  used. 

The  plotting  space  available  with  DRAW  is  fifteen  inches  in 
the  y-direction  and  nine  inches  in  the  x-direction.   Recall  that 
the  convention  for  numbering  the  nodes  usually  results  in  having 
the  narrowest  side  of  the  structure  along  the  y-axis  (section  3. A). 
Therefore,  GRID  interchanges  the  x-  and  y-coordinates  of  the  structure, 
for  plotting  purposes  only,  to  take  advantage  of  the  larger  plotting 
space  in  the  y-direction.   This  further  reduces  the  scale  of  the 
plot. 
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IV.   INPUT  DATA  PREPARATION 

This  chapter  is  written  as  a  self-contained  unit  giving  specific 
information  necessary  for  the  preparation  of  input  data  cards.   The 
chapter  may  be  removed  and  used  as  a  guideline  for  input  data  preparation 
after  a  thorough  knowledge  of  the  computer  programs  has  been  obtained. 

A.   VARIABLES 

The  variables  and  formulas  defined  below  are  included  in  the 
code  and  should  be  used  as  an  aid  for  determining  the  arrangement 
of  nodes  within  the  mesh  when  preparing  input  data. 
Variable  Meaning 

1.  NEL  Total  number  of  elements. 

NEL  =  (NX) (NY) (NZ)  ,  where  NX,  NY  and  NZ  are  the  number  of  elements 
along  the  §  ,T)  ,  and  Q     directions     (21) 

2.  NPS  Total  number  of  nodes  on  any  face  at  constant 

£  ,  where  a  face  is  a  slice  which  intersects 
both  corner  nodes  and  mid-side  nodes. 

NPS  -  3 (NX) (NY)  +  2(NX  +  NY)  +  1      (QUAMEG)  (22) 

NPS  =  5 (NX) (NY)  +  3 (NX  +  NY)  +1      (CUMEG)  -   (23) 

3.  NPM  Total  number  of  mid-side  nodes  between  any  two 

faces. 

NPM  -  (NX  +  1) (NY  +  1)        (QUAMEG)  (24) 

NPM  =  2  (NX  +  1)(NY  +  1)       (CUMEG)  (25) 

4.  NPL  Total  number  of  nodes  in  one  layer,  where  a 

layer  consists  of  one  face  plus  the  slice 
intersecting  mid-side  nodes  between  that  face 
and  an  adjacent  face. 

NPL  =  NPM  +  NPS  (26) 

5.  NUMNP  Total  number  of  nodes  in  the  mesh. 

NUMNP  -  (NZ)  (NPL)  +  NPS  (27) 
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B>   DATA  CARDS 

L   El^enL^i^lOa^a  (7110,215) 

Number  of  cards  per  problem  -  one 


Columns 

Variable 

1  -  10 

NX 

11  -  20 

NY 

21  -  30 

NZ 

31  -  40 

NDX 

41  -  50 

NDY 

51  -  60 

NDZ 

61  -  70 

NPDP 

Meaning 


71  -  75 


KCRB 


Number  o£  elements  in  the  ^'-direction; 
Number  of  elements  in  the  H'-Olrection; 
Number  of  elements  in  the  C'-ilrections 
Number  of  segments  in  the  ^'-direction; 
Number  of  segments  in  the  V-direction; 
Number  of  segments  in  the  c'-dlrection; 

„f  nodes  for  which  the  coor- 
Total  number  of  nodes  *      ^.^  ^ 

dinates  must  be  spetu.  other 

location  of  segment  end  poxnts  or 
boundary  nodes ; 
s   ifieswhetber  rectangular^or^lindri. 

Ca  L     LWiui  i„Ai  rates     CV  LinuL  !.«-«- ~ 

A  positive  i^f  \  f^Cao  e  e  I  specifies 
coordinates  and  a  blank  or  zero  v 
rectangular  coordinates; 

2.   rrnhl-TTi  Identification  (10A8) 

Number  of  cards  per  problem  -  one 

Meaning 
Columns      Variable 

-—      TITLE       Alphameric  information  to  identrfy  the 
*  "  bU  problem. 

3.   No^^oint^c^rdinates.  (1110,  3F10.0) 

..      onP  for  each  NPDP  point 
Number  of  cards  per  problem  -  one 

Meaning 
Columns      Variable 

Nodal  point  number; 
1-10      1 
u.20      COORD(I.l)    x-coordinate  (radius); 


/9 


Columns 


Variable 


Meaning 
y-coordinate  (angle) ; 
z-coordinate . 


21  -  30      C00RD(I,2) 
31  -  40      COORD (I, 3) 

4.   Segment  Divisions   (5110) 
Number  of  cards  per  problem  -  one  for  each  segment 


Columns 

Variable 

1-10 

I 

11  -  20 

J 

21  -  30 

Nl 

31  -  40 


Ml 


41  -  50 


LI 


Meaning 

Lowest  node  number  in  the  segment; 

Highest  node  number  in  the  segment; 

Segments  in  the  §'-  and  ^'-directions: 
Number  of  elements  above  the  segment 
(in  the  positive  Tj'-direction)  ;  Segments 
in  the  T|'-direction:   Number  of  elements 
within  the  segment; 

Segments  in  the  T|'-  and  Q   -directions: 
Number  of  elements  to  the  left  of  the 
segment  (in  the  negative  §  -direction); 
Segments  in  the  §  -diie<-  Liori ;   Number 
of  elements  within  the  segment; 

Segments  in  the  ?'-  and  T]'-directions  : 
Number  of  elements  before  the  segment 
(in  the  positive  G'-direction) ;   Segments 
in  the  (/-direction:   Number  of  elements 
within  the  segment. 


5.  Load  Data  (2110,  2F10.0) 
Number  of  cards  per  problem  -  one 
Columns      Variable 


Meaning 


1-10 


MAP 


11  -  20 
21  -  30 
31  -  40 


NEFP 

SGZ 

UDP 


A  positive  integer  indicates  that  a 
plot  of  the  structure  mesh  is  desired. 
A  blank  or  zero  indicates  that  a  plot 
is  not  desired; 

Number  of  pressurized  element  faces; 

Specific  weight ; 

Uniformly  distributed  pressure. 
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6.   Pressure  Loading  Data  (2110) 

Number  of  cards  per  problem  -  one  for  each  NEFP 

Columns       Variable  Meaning. 

1-10       NELP  Element  number  on  which  a  pressure  exists; 

u  _  20       NFACE  Face  on  which  the  pressure  is  applied. 

(See  section  3.G  and  Figure  6) 

7 .   Plot  Title  Cards   (6A8) 
Number  of  cards  per  problem  -  two  per  plot  (NZ  +  1  plots) 
Columns      Variable  Meaning. 

1  _  48       ITITLE  Alphameric  information  to  identify  the 

plot; 

,  _  43       ITITLE  Alphameric  information  to  identify  the 

programmer  and  his  computer  center  box 
number . 

8.   Stop  Trap 

The  stop  trap  employed  in  the  programs  allows  the  execution 
of  as  many  separate  problems  as  may  be  desired  with  one  submission 
of  the  program.   The  program  will  continue  processing  problems  until 
it  is  informed  that  NX  is  less  than  zero.   Therefore,  the  last  card 
of  a  complete  input  data  deck  should  have  a  negative  integer  in 
columns  1-10. 
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V.   CONCLUSIONS 

The  addition  of  QUAMEG  and  CUMEG  to  the  tools  available  in  the 
finite  element  method  of  analysis  opens  new  doors  for  the  analyst 
who  uses  isoparametric  finite  elements.   Previously,  the  amount 
of  input  data  required  to  analyze  a  structure  of  appreciable  size 
was  too  overwhelming  for  the  analyst  to  undertake  the  problem. 
QUAMEG  and  CUMEG  make  it  practical  for  the  analyst  to  solve  such 
a  problem.   In  addition,  he  can  be  certain,  after  an  acceptable 
mesh  is  generated,  that  the  output  data  generated  by  QUAMEG  and 
CUMEG  is  free  of  error. 

For  the  example  discussed  in  Appendix  A,  the  number  of  input 
cards  required  by  QUAMEG  versus  the  number  of  output  cards  generated 
by  QUAMEG  tor  direct  input  to  TRISOP  is  about  six  per  cent.   For 
the  1  x  8  x  8  mesh  of  the  pinched  cylinder  problem  discussed  in 
Appendix  B,  this  ratio  is  about  4  per  cent. 

With  a  gravity  load  and  pressure  load  acting  on  one  cubic 
element,  it  was  found  that  QUAMEG  and  CUMEG  produced  consistent 
load  vectors  which  were  accurate  to  at  least  eight  decimal  places. 

QUAMEG  and  CUMEG  considerably  reduce  the  effort  required  of 
the  analyst  to  solve  a  problem  by  the  finite  element  method,  generate 
acceptable  values  for  nodal  point  coordinates  and  compute  accurate 
consistent  gravity  and  pressure  load  vectors. 
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VI.      RECOMMENDATIONS 

A.  LOADING   CONDITIONS 

The  capabilities  of  QUAMEG  and  CUMEG  should  be  increased  by 
adding  subroutines  which  will  generate  the  consistent  load  vectors 
for  a  thermal  gradient  in  the  structure  and  for  an  initial  displace- 
ment condition  at  one  or  more  nodes. 

B.  PLOT  OF  STRUCTURE  MESH 

It  would  be  preferred  to  obtain  a  more  complete  representation 
of  the  structure  mesh  than  is  now  provided.  Two  courses  of  action 
are  recommended : 

1.  Expand  Subroutine  GRID  in  such  a  way  that  plots  of  faces 
in  the  T)  -C,     and/or  §  -Q'    planes  are  also  generated. 

2.  Replace  GRID  with  a  subroutine  which  would  generate  a 
three-dimensional  plot  of  the  structure  mesh.   The  feasibility 
of  using  the  IBM  360  source  library  program  PLOT  PACKAGE,  which 
is  mentioned  in  Ref.  12,  should  be  investigated  if  this  method 
is  attempted. 

This  author  feels  that  the  results  of  recommendation  2  would 
be  preferred  to  1,  although  the  coding  may  be  more  difficult. 
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APPENDIX  A 
GENERATING  A  MESH  WITH  Q1JAMEG 


This  appendix  describes  the  preparation  of  input  cards  and  the 
expected  results  for  a  sample  problem  using  QUAMEG.   This  problem 
was  chosen  as  an  example  to  demonstrate  the  capabilities  of  QUAMEG 
and  to  more  fully  explain  the  procedures  described  in  Chapter  4. 

The  structure  considered  is  a  cylindrical  solid  with  a  star 
shaped  void  in  the  center.   This  void  linearly  converges  to  a  point. 
Figure  7  shows  a  45   sector  of  the  structure.   It  is  normally  prefer- 
able to  work  with  angles  in  the  first  quadrant.   Therefore,  the 
sector  from  zero . to  forty- five  degrees  will  be  considered. 

The  first  step  in  ulie  problem  is  to  decide  upon  a  mesh  size. 
In  this  case,  an  8  x  5  x  4  mesh  was  chosen.   The  amount  of  control 
to  be  exercised  by  the  programmer  to  insure  that  the  structure 
will  be  modeled  accurately  can  be  determined  now.   The  coordinates 
of  all  nodes  along  the  edges  of  the  bounding  surfaces  of  the  structure 
must  be  available  to  the  system  before  the  remaining  nodal  coor- 
dinates can  be  determined.   Since  there  are  no  irregularities  on 
the  boundaries  between  the  front  and  back  faces ,  the  coordinates 
of  the  nodes  on  the  edges  between  these  two  faces  can  be  computed 
by  providing  the  necessary  segment  data.   Therefore,  only  the  nodes 
on  the  front  and  back  faces  need  be  considered. 

At  this  point,  two  rectangles  of  constant  Q      are  drawn  with 
the  proposed  mesh  superimposed  upon  them  (see  Figures  8  and  9) . 
These  figures  will  serve  as  an  aid  for  numbering  the  nodes  on  the 
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actual  structure  and  for  determining  the  locations  of  segments. 
After  the  nodes  on  face  one  are  numbered,  equation  22  is  solved  to 
insure  that  an  error  has  not  been  made  in  the  numbering.   The  result 
of  equation  22  is  subtracted  from  that  of  equation  27  to  determine 
the  number  of  the  first  node  on  the  back  face.   In  this  way,  the 
programmer  need  not  go  through  the  time-consuming  process  of  numbering 
all  nodes  in  the  mesh.   The  nodes  on  the  back  face  are  numbered  and 
checked  as  before.   Note  that  the  8  x  5  mesh  is  in  the  S'-T)'  plane. 
Orienting  the  structure  so  that  the  5  x  4  face  of  the  mesh  were  in 
the  g'-V  plane  would  reduce  the  half-bandwidth  of  the  system. 
However,  cylindrical  coordinates  must  be  used  in  this  problem  and 
it  would  not  be  possible  to  model  the  structure  accurately  if  it  were 
oriented  in  that  manner.   Finally,  a  picture  is  drawn  of  the  front 
and  back  faces  of  the  structure  with  the  desired  mesh  superimposed 
and  the  nodes  numbered  (Figures  10  and  11) . 

The  identification  of  segments  will  now  be  determined.   The 
segments  will  be  considered  in  the  same  order  as  they  are  considered 
in  the  program.   Referring  to  Figures  8  and  9,  it  is  seen  that  the 
following  lines  are  parallel  to  the  Tl'-axis:   1  -  11,  137  -  147, 
805  -  815,  941  -  951.   Now  referring  to  Figures  10  and  11,  it  is 
seen  that  all  of  these  lines  are  of  constant  angle  with  a  uniform 
distribution  of  elements  within  the  lines.   Therefore,  these  lines 
will  be  treated  as  segments  and  the  coordinates  of  the  nodes  listed 
above  will  be  included  in  the  input  data. 

By  this  method  of  specifying  lines,  the  author  does  not  intend 
to  imply  that  the  nodes  along  the  line  are  numbered  consecutively, 
starting  with  the  lower  number.   For  example,  the  nodes  along  line 
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60  -  111  in  Figure  8  are:   60,  67,  77,  84,  94,  101,  and  111.   This 
method  of  specifying  lines  merely  lists  the  numbers  of  nodes  at  the 
ends  of  a  line. 

In  order  to  control  the  departure  angles  of  lines  79  -  77, 
96  -  94,  and  113  -  111,  it  is  necessary  to  force  line  52  -  62  nearer 
the  x-axis  than  would  result  if  an  even  distribution  of  nodes  along 
line  1  -  137  were  allowed.   Therefore,  an  appropriate  location  for 
node  52  is  chosen,  the  coordinates  of  nodes  52  and  62  are  included 
in  the  input  data,  and  line  52  -  62  is  treated  as  a  segment.   Note 
that  this  line  is  of  neither  constant  angle  nor  constant  radius. 
However,  since  the  line  is  not  on  the  edge  of  a  bounding  surface, 
the  fact  that  it  will  be  curved  will  not  affect  the  solution  adversely 
Similarly,  the  coordinates  of  nodes  856  and  866  will  be  included 
in  the  input  data  and  that  line  will  be  a  segment. 

The  same  rules  govern  the  selection  of  segments  in  the  |  -direc- 
tion.  Referring  to  Figures  8  and  9,  it  is  seen  that  the  coordinates 
of  nodes  along  lines  1  -  137 ,  11  -  147,  805  -  941.  and  815  -  951 
must  be  made  available  to  the  computer.   In  Figures  10  and  11,  line 

I  -  137  is  a  line  of  constant  radius;  however,  the  element  divisions 
along  this  line  are  not  uniform.   Two  segments  must  be  specified: 
lines  1-52  and  52  -  137.   The  coordinates  of  the  end  nodes  on  these 
segments  have  been  specified  already  and  need  not  be  repeated.   Line 

II  -  62  is  of  constant  radius  and  is  therefore  one  segment.   Again, 
the  coordinates  of  these  nodes  have  been  specified  previously.   Line 
62  -  113  is  of  constant  angle,  meeting  the  criterion  for  a  segment. 
The  coordinates  of  node  113  must  be  specified.   Line  113  -  147  is 

of  neither  constant  angle  nor  constant  radius.   Therefore,  it  cannot 
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Pnt  and  the  coordinates  of  all  nodes  along  this 
be  treated  as  a  segment  and  tne  co 

-  A»+a        Tn  order  to  further  control 
line  „st  be  Included  in  the  input  data. 

the  departure  angle  of  the  lines  mentioned  ahove,  line  60  -  145 
ls  also  included  as  a  segment.   The  coordinates  of  the  end  nodes 

roH  in  the  problem  already  and  do 
,  .       ~r,f  havp  been  computed  in  tne  viu 
of  this  segment  nave  uecu     v 

not  need  to  he  included  in  the  input  data. 

the  exact  shape  of  line  60  -  W  is  not  important.   Similarly  on 

~,  -i     a    QAA  -  949  are  treated 
the  bac.  face,  lines  805  -  856,  856  -  941,  and  864 

ents   Although  the  nodes  along  line  815  -  95!  are  at  the 
as  segments.   rtLLUUU& 

origin,  care  must  he  ta.en  to  insure  that  the .coordinates  of  these 

nodes  are  specified  at  appropriate  angles  in  such  a  way  that  this 

m  us,    p"  as  the  front  face.   The  radial 
point  will  have  the  same  "shape   as 

a-    afPs  is  zero.   This  care  is  necessary 
component  of  these  coordinates  is  zero. 

because  all  computations  in  the  problem  are  made  with  cylindrical 

coordinates.  Were  the  coordinates  of  the  nodes  along  this  line 

„lo  nf  7ero   the  mesh  would 
specified  at  a  radius  of  zero  and  an  angle  of  zero, 

be  distorted.   The  following  lines  can  also  he  treated  as  segments: 
815-  866  and  866  -917.   The  coordinates  of  nodes  923,  934,  and 

u   •    i-  Hsra  for  the  same  reason  as  was 
940  must  be  included  in  the  input  data 

given  for  nodes  119,  130,  and  136. 

To  complete  the  process  of  specify  the  coordinates  of  nodes 
along  the  structure  edges,  the  following  lines  will  be  treated  as 

ents  in  the  Cation..   1-805,52-856   137-941,60-864, 
segments  m  me  s, 
77.  331,  94-  898,  11  -915,  11  -815,  62  -866,  113  -917,  and 


130  -  934. 

The  procedure  for  determining  the  values  of  *l.  Ml.  and  LI 
for  the  segment  division  data  will  not  be  discussed  here.   It  is 
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felt  that  the  explanation  in  Chapter  4,  combined  with  the  listing 
of  the  input  cards  for  this  problem,  is  sufficiently  clear. 

The  card  containing  the  element  field  data  is  prepared  after 
the  above  steps  have  been  accomplished.   Loading  information  and  the 
request  for  punched  cards  should  not  be  included  in  the  input  data 
until  it  has  been  determined  that  the  mesh  will  be  satisfactory. 

Approximately  45  seconds  of  computer  time  were  required  to  solve 
this  problem  with  QUAMEG,   That  time  includes  the  31  seconds  required 
to  compile  the  program.   Figures  12,  13,  14,  15,  and  16  are  the  plots 
of  the  structure  mesh  that  was  generated  by  QUAMEG .   Table  I  is  the 
listing  of  the  required  input  cards  described  above.   Figure  17  is 
the  result  of  a  5  x  5  x  3  mesh  generated  with  CUMEG  for  the  same 
structure.   Intermediate  slices  and  the  back  surface  are  not  shown. 
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Table  I 
Listing  of  Input  Data  Cards 
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Figure  12.   Generated  Mesh  For  Front  Face 
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Figure  13.   Generated  Mesh  For  First  Intermediate  Slice 
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Figure  14.   Generated  Mesh  For  Second  Intermediate  Slice 
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Figure  15.   Generated  Mesh  For  Third  Intermediate  Slice 
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Figure  16.   Generated  Mesh  For  Back  Face 
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Figure  17.   Generated  Mesh  For  Front  Face 
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APPENDIX  B 
A  CLASSICAL  PROBLEM  SOLVED  WITH  TRISOP 


TRISOP  is  a  relatively  new  computer  system  which  uses  isoparametric 
finite  elements  to  obtain  stresses,  strains,  and  displacements  of 
three-dimensional  structures  that  are  subjected  to  static  loads. 

Although  a  number  of  test  problems  have  been  conducted  with 
TRISOP,  the  input  preparation  was  too  voluminous  to  contemplate 
a  full  evaluation  of  the  system.   Without  QUAMEG ,  the  difficulties 
discussed  in  Chapter  1  limited  the  size  and  number  of  problems  that 
could  be  solved.   In  order  to  evaluate  TRISOP  as  a  complete  system, 
a  classical  problem  was  selected  and  a  convergence  study  was  made, 
using  QUAMEG  and  TRISOP. 

The  problem  selected  was  that  of  the  pinched,  thin-shell  cylinder. 
Obviously,  this  problem  should  be  solved  using  thin-shell  elements 
for  maximum  efficiency.   The  purpose  of  this  test  was  to  determine 
if  the  three-dimensional  solution  would  converge  to  the  same  results 
as  are  reported  in  Ref.  3  and  also  to  determine  how  well  thin-shell 
action  is  captured  by  TRISOP.   The  results  reported  by  Cantin  [3] 
are  included  in  Table  II  with  the  results  obtained  in- the  test.   The 
problem  is  described  in  Figure  18. 

Conclusions : 

Although  the  results  of  this  convergence  study  indicate  a 
tendency  to  approach  the  same  results  obtained  by  thin-shell  theory, 
acceptable  results  can  be  obtained  only  with  a  very  fine  mesh. 
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TRISOP  can  be  used  to  solve  problems  for  which  there  are  no 
closed-form  or  accurate  numerical  methods  of  solution  available. 
However,  the  user  must  be  prepared  to  pay  the  price  of  long  compu- 
tational times,  for  some  problems. 


Table  II 
Displacement  at  the  applied  load  for  a  pinched  cylinder. 


Bogner  et  al. 


Cantin 


No. 

of 

Mesh 

Displace. , 

Eq. 

in. 
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1x2 

-0.0802 

96 

1x3 

-0.1026 
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1  x  4 

-0.1087 

108 

2x2 
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2x3 

-0.1036 

180 

2  x  k 

-0.1098 

No. 

of 

Mesh 

Displac . , 

Eq. 

in. 

36 

1x2 

-0.0921 

48 

1  x  3 

-0.1072 

60 

1x4 

-0.1099 

54 

2x2 

-0.0931 

7  2 

2x3 
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90 

2x4 

-0.1113 
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4x4 

-0.1126 
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6x6 

-0.1137 
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8x8 

-0.1139 

7  26 


10  x  10   -0.1139 


Hanson 

No. 

of 

Eq. 
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Displac . , 
in . 

153 

243 
423 
465 

945 

1 
2 
2 

1 
1 

x  2  x  2 
x  2  x  2 
x  2  x  4 
x  4  x  4 
x  6  x  6 

-0.00177 
-0.00178 
-0.00196 
-0.00375 
-0.03002 

1593 


1x8x8 


-0.08702 
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APPEND IX  C 
COMPUTER  LISTING  (QUAMEG) 
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APPENDIX  D 
COMPUTER  LISTING  (CUMEG) 
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